# Replication code for "Loss Framing in Territorial Disputes"
# Andi Zhou, H.E. Goemans, and Michael Weintraub
# January 2025
# R version 4.2.2 (2022-10-31)
# RStudio version 2024.04.1+748
# See readme for version numbers of all packages used

library(tidyverse)
library(stargazer)
library(ggplot2)
library(mediation)
library(margins)
library(cobalt)
library(ivreg)
library(interflex)
library(car)

# Set to working directory of choice
# setwd()

# Load datasets
arg_dat <- read.csv("arg_rep.csv")
chl_dat <- read.csv("chl_rep.csv")


#################
### MAIN TEXT ###
#################

### Fig. 4, Panel A: Argentina main results ###
# Prepare results to plot
arg_summary <- arg_dat %>% 
  mutate(opp_name = factor(opp_name, c("No opponent", "Chile", "UK"))) %>%
  group_by(opp_name, frame_loss) %>% 
  summarise(mean_gamble = mean(gamble, na.rm = T),
            se = sd(gamble, na.rm = T)/sqrt(n() - 1),
            lo_CI = mean_gamble - 1.96*se,
            hi_CI = mean_gamble + 1.96*se)

# Create plot
ggplot(data = arg_summary, aes(x = opp_name, y = mean_gamble, color = as.factor(1-frame_loss))) +
  geom_point(size = 2)+
  geom_text(aes(label = round(mean_gamble, digits = 3)), nudge_x = 0.25, show.legend = F)+
  geom_errorbar(aes(ymin = lo_CI, ymax = hi_CI), width = 0.1)+
  theme(plot.title = element_text(hjust = 0.5))+
  ylim(c(0, 1))+
  ylab("Proportion choosing risky policy option")+
  xlab("Opponent")+
  theme(legend.position = "bottom")+
  scale_color_manual(values = c("gray45", "black"), labels = c("Loss", "Gain"))+
  labs(color = "")

### Fig. 4, Panel B: Chile main results ###
# Prepare results to plot
chl_summary_opp <- chl_dat %>% 
  drop_na(opp_name) %>%
  group_by(opp_name, frame_loss) %>% 
  summarise(mean_gamble = mean(gamble, na.rm = T),
            se = sd(gamble, na.rm = T)/sqrt(n() - 1),
            lo_CI = mean_gamble - 1.96*se,
            hi_CI = mean_gamble + 1.96*se)

# Create plot
ggplot(data = chl_summary_opp, aes(x = opp_name, color = as.factor(1-frame_loss))) +
  geom_point(aes(y = mean_gamble), size = 2)+
  geom_text(aes(y = mean_gamble, label = round(mean_gamble, digits = 3)), nudge_x = 0.2, show.legend = F)+
  geom_errorbar(aes(ymin = lo_CI, ymax = hi_CI), width = 0.1)+
  theme(plot.title = element_text(hjust = 0.5))+
  scale_color_manual(values = c("gray45", "black"), labels = c("Loss", "Gain"))+
  ylim(c(0, 1))+
  ylab("Proportion choosing risky policy option")+
  xlab("Opponent")+
  theme(legend.position = "bottom")+
  labs(color = "")

## Footnote 15: Linear hypothesis test showing no significant difference between Chile vs. no opponent effect and UK vs. Chile effect
# Estimate OLS model, interacting loss framing treatment with opponent treatments
opp_effect_arg <- lm(gamble ~ frame_loss + any_opp + uk + frame_loss:any_opp + frame_loss:uk, arg_dat)

# Check whether the uninteracted opponent effect coefficients (i.e. the opponent effects in gain frame) are equal to each other
linearHypothesis(opp_effect_arg, "any_opp = uk")


### Fig. 5: Main results ###
# Prepare results to plot
summary_joint <- chl_dat %>% drop_na(historical) %>%
  group_by(hi_val, historical, frame_loss) %>% 
  summarise(mean_gamble = mean(gamble, na.rm = T),
            se = sd(gamble, na.rm = T)/sqrt(n() - 1),
            lo_CI = mean_gamble - 1.96*se,
            hi_CI = mean_gamble + 1.96*se) %>% 
  mutate(condition = paste0(c("Low val.", "High val.")[hi_val + 1], ", ", c("Non-hist.", "Historical")[historical + 1]))

# Create plot
ggplot(data = summary_joint, aes(x = condition, color = as.factor(1-frame_loss))) +
  geom_point(aes(y = mean_gamble), size = 2)+
  geom_text(aes(y = mean_gamble, label = round(mean_gamble, digits = 3)), nudge_x = 0.2, show.legend = F)+
  geom_errorbar(aes(ymin = lo_CI, ymax = hi_CI), width = 0.1)+
  theme(plot.title = element_text(hjust = 0.5))+
  scale_color_manual(values = c("gray45", "black"), labels = c("Loss", "Gain"))+
  ylim(c(0, 1))+
  ylab("Proportion choosing risky policy option")+
  xlab("")+
  theme(legend.position = "bottom")+
  guides(color = guide_legend(reverse=TRUE)) +
  labs(color = "")+
  coord_flip()


### Fig. 6: Marginal effect of historical ownership treatment by framing and economic value conditions ###
## Estimate logit model with framing, historical ownership, and economic value treatments fully interacted plus demographic and dispositional controls
mod_hist_logit_marg <- glm(gamble ~ frame_loss*historical*hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "logit"), data = chl_dat)

## Generate marginal effects from logit model
# Gain frame
hist_marg_val_gain <- cplot(mod_hist_logit_marg, x = "hi_val", dx = "historical", what = "effect", 
                            data = chl_dat %>% filter(frame_loss == 0), draw = F, n = 2) %>% data.frame()
# Loss frame
hist_marg_val_loss <- cplot(mod_hist_logit_marg, x = "hi_val", dx = "historical", what = "effect", 
                            data = chl_dat %>% filter(frame_loss == 1), draw = F, n = 2) %>% data.frame()

## Prepare results to plot
hist_marg_val_gain$frame <- "Gain frame"
hist_marg_val_loss$frame <- "Loss frame"
hist_marg_val_gain$condition <- hist_marg_val_loss$condition <- factor(c("Low value", "High value"), 
                                                                       levels = c("Low value", "High value"))
## Bind gain-frame and loss-frame results
hist_marg_val <- rbind(hist_marg_val_gain, hist_marg_val_loss)

## Create plot
ggplot(hist_marg_val, aes(x = condition, y = yvals))+
  geom_pointrange(aes(ymin = lower, ymax = upper)) +
  geom_text(aes(label = round(yvals, digits = 3)), nudge_x = 0.25)+
  geom_hline(yintercept = 0, linetype = "dashed")+
  theme(plot.title = element_text(hjust = 0.5))+
  xlab("")+
  ylab("Marginal effect of historical framing")+
  facet_wrap(~frame)+
  theme(legend.position = "none")


### Z-score of most significant effect (gain frame, high value condition): 1.766, or p = 0.077 (cited on p. 21) ###
hist_marg_val$yvals[1]/((hist_marg_val$upper[1] - hist_marg_val$yvals[1])/1.959)


### Likelihood ratio test shows historical, value variables add no explanatory power (cited on p. 21) ###

# Define parsimonious model without historical ownership and economic value treatments
mod_hist_logit_pars <- glm(gamble ~ frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "logit"), data = subset(chl_dat, !is.na(historical) & !is.na(hi_val)))

# Conduct likelihood ratio test against the model used for Figure 6
anova(mod_hist_logit_pars, mod_hist_logit_marg, test = "Chisq")


### Fig. 7: Results with perceived ownership DV ###

## Left-hand panel ##
chl_dat %>% filter(!is.na(historical), !is.na(perceive_own_lab)) %>% 
  group_by(historical, hi_val, perceive_own_lab) %>% 
  summarise(value = n()) %>%
  mutate(condition = paste0(c("Non-historical", "Historical")[historical+1], ", ",
                            c("low value", "high value")[hi_val+1])) %>%
  ggplot(aes(fill = factor(perceive_own_lab, levels = c("Not sure", "No", "Yes")), y = value, x = condition)) +
  geom_bar(position = "fill", stat = "identity", width = 0.5) +
  theme_bw()+
  labs(fill ="")+
  ylab("% perceive Chile owns")+
  xlab("Treatments")+
  theme(legend.position = "bottom")+
  guides(fill = guide_legend(reverse=TRUE)) +
  scale_fill_manual(values = c("gray", "gray45", "black"),
                    labels = c("Not sure", "No", "Yes"))+
  coord_flip()


## Right-hand panel ##
# Difference of means, low value condition
diff_lo_val <- mean(subset(chl_dat, historical == 1 & hi_val == 0)$perceive_own_bin, na.rm = T) - 
  mean(subset(chl_dat, historical == 0 & hi_val == 0)$perceive_own_bin, na.rm = T)
# Difference of means, high value condition
diff_hi_val <- mean(subset(chl_dat, historical == 1 & hi_val == 1)$perceive_own_bin, na.rm = T) - 
  mean(subset(chl_dat, historical == 0 & hi_val == 1)$perceive_own_bin, na.rm = T)

# T-test, low value condition
t_lo_val <- t.test(subset(chl_dat, historical == 1 & hi_val == 0)$perceive_own_bin, subset(chl_dat, historical == 0 & hi_val == 0)$perceive_own_bin)
# T-test, high value condition
t_hi_val <- t.test(subset(chl_dat, historical == 1 & hi_val == 1)$perceive_own_bin, subset(chl_dat, historical == 0 & hi_val == 1)$perceive_own_bin)

# Combine results in table for plotting
perceive_own_plot <- data.frame(rbind(c(diff_lo_val, t_lo_val$conf.int),
                           c(diff_hi_val, t_lo_val$conf.int)))

# Add column names
colnames(perceive_own_plot) <- c("diff", "lo_CI", "hi_CI")

# Add labels for economic value conditions
perceive_own_plot$condition <- c("Low value", "High value")

# Create plot
ggplot(perceive_own_plot, aes(x = condition))+
  geom_pointrange(aes(y = diff, ymin = lo_CI, ymax = hi_CI))+
  geom_text(aes(y = diff, label = round(diff, digits = 4)), nudge_x = 0.1)+
  ylim(c(-.05, .2))+
  geom_hline(aes(yintercept = 0), linetype = "dashed")+
  ylab("Diff. in means (Historical - non-historical)")+
  xlab("Value treatment")+
  coord_flip()


### Table 2 ###
list(lm(force_ab ~ gamble + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     ivreg(force_ab ~ gamble + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
             frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     ivreg(force_ab ~ gamble + policy_mil + I(gamble * policy_mil) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
             frame_loss + policy_mil + I(frame_loss * policy_mil) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     ivreg(force_ab ~ gamble + policy_mil + I(gamble * policy_mil) + I(gamble * mi_scale) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
             frame_loss + policy_mil + I(frame_loss * policy_mil) + I(frame_loss * mi_scale) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat)) %>%
stargazer(type = "text",
          table.placement = "!htbp",
          title = "Risk acceptance on support for military action, OLS and IV results. Models 2 through 4 use loss framing as the instrumental variable for risk acceptance.",
          label = "tab:force",
          omit.stat = c("adj.rsq", "ser", "f"),
          dep.var.labels = "Support military action",
          column.labels = c(rep("\\textit{OLS}", 1), rep("\\textit{IV}", 3)),
          model.names = F,
          order = c("gamble", "gamble*", "policy", "mi_scale"),
          covariate.labels=c("Risky option", "Risky option x Mil. action", "Risky option x Mil. assert.", "Mil. action", "Mil. assertiveness", "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", "National identif.", "Int'l trust", "Social trust"),
          no.space = T)



#########################
### ONLINE APPENDICES ###
#########################

### Table A1 ###
## Create table
arg_bal_tab_opp <- arg_dat %>%
  # Group by treatment cells
  group_by(frame_loss, opp_name) %>%
  # Get mean value of each demographic variable by treatment cell
  summarise(across(c(age, female, race_white, university, buenos_aires, ideology_scale), ~ mean(.x, na.rm = T))) %>%
  ungroup() %>% dplyr::select(age:ideology_scale) %>% pivot_longer(cols = age:ideology_scale) %>% group_by(name) %>%
  # Get minimum and maximum mean values of each demographic variable
  summarise(min = min(value), max = max(value)) %>%
  # Re-order table
  mutate(name = factor(name, levels = c("age", "female", "race_white", "university", "buenos_aires", "ideology_scale"))) %>% arrange(name) %>%
  # Merge population means of each demographic variable
  mutate(target_means = c(41.84, # Age; LAPOP 2018
                          .504, # Female; LAPOP 2018
                          .540, # White; LAPOP 2018
                          0.337, # University education; LAPOP 2018
                          .462, # Buenos Aires residents as share of population; 2010 Census
                          5.204), # Ideology; LAPOP 2018
         # Merge overall sample mean of each demographic variable
         sample_means = arg_dat %>% dplyr::select(age, female, race_white, university, buenos_aires, ideology_scale) %>% colMeans(na.rm = T))

## Calculate standardized mean differences using 'cobalt' package
# Drop all observations with missing demographic data to use in cobalt
arg_dat_sub <- arg_dat %>% dplyr::select(frame_loss, opp_name, gamble, age, female, race_white, university, buenos_aires, ideology_scale) %>% na.omit

# Estimate standardized mean differences using cobalt
arg_cobalt <- bal.tab(arg_dat_sub %>% dplyr::select(age:ideology_scale), treat = paste0(arg_dat_sub$frame_loss, arg_dat_sub$opp_name))

# Incorporate standardized mean difference estimates from cobalt into table
arg_bal_tab_opp$std_mean_diff <- arg_cobalt$Balance.Across.Pairs[,2] %>% round(3)

## Clean up 
# Reorder columns
arg_bal_tab_opp <- dplyr::select(arg_bal_tab_opp, name, target_means, sample_means, max, min, std_mean_diff)

# Label rows and columns
arg_bal_tab_opp$name <- c("Age", "Female", "White", "University", "Buenos Aires", "Right-wing ideology (1-10)")
colnames(arg_bal_tab_opp) <- c("Covariate", "Pop. mean", "Sample mean", "Max. cell mean", "Min. cell mean", "Std. mean diff.")

## Print table
print(arg_bal_tab_opp)


### Table A2 ###
## Create table
chl_bal_tab_opp <- chl_dat %>%
  # Group by treatment cells
  group_by(frame_loss, opp_uk) %>%
  # Get mean value of each demographic variable by treatment cell
  summarise(across(c(age, female, race_white, educ_superior, santiago, ideology_scale), ~ mean(.x, na.rm = T))) %>%
  ungroup() %>% dplyr::select(age:ideology_scale) %>% pivot_longer(cols = age:ideology_scale) %>% group_by(name) %>%
  # Get minimum and maximum mean values of each demographic variable
  summarise(min = min(value), max = max(value)) %>%
  # Re-order table
  mutate(name = factor(name, levels = c("age", "female", "race_white", "educ_superior", "santiago", "ideology_scale"))) %>% arrange(name) %>%
  # Merge population means of each demographic variable
  mutate(target_means = c(42.38, # Age; LAPOP 2021 
                          8946991/17520239, # Female; Census 2017
                          0.437, # White; LAPOP 2021
                          0.298, # Superior educ.; Census 2017
                          7112808/17520239, # Santiago residents as share of population; Census 2017
                          5.054), # Ideology; LAPOP 2018
         # Merge overall sample mean of each demographic variable
         sample_means = chl_dat %>% filter(!is.na(opp_uk)) %>% dplyr::select(age, female, race_white, educ_superior, santiago, ideology_scale) %>% colMeans(na.rm = T))

## Calculate standardized mean differences using 'cobalt' package
# Drop all observations with missing demographic data to use in cobalt
chl_dat_sub <- chl_dat %>% dplyr::select(frame_loss, opp_uk, gamble, age, female, race_white, educ_superior, santiago, ideology_scale) %>% na.omit

# Estimate standardized mean differences using cobalt
chl_cobalt <- bal.tab(chl_dat_sub %>% dplyr::select(age:ideology_scale), treat = paste0(chl_dat_sub$frame_loss, chl_dat_sub$opp_uk))

# Incorporate standardized mean difference estimates from cobalt into table
chl_bal_tab_opp$std_mean_diff <- chl_cobalt$Balance.Across.Pairs[,2] %>% round(3)

## Clean up 
# Reorder columns
chl_bal_tab_opp <- dplyr::select(chl_bal_tab_opp, name, target_means, sample_means, max, min, std_mean_diff)

# Label rows and columns
chl_bal_tab_opp$name <- c("Age", "Female", "White", "Superior ed.", "Santiago", "Right-wing ideology (1-10)")
colnames(chl_bal_tab_opp) <- c("Covariate", "Pop. mean", "Sample mean", "Max. cell mean", "Min. cell mean", "Std. mean diff.")

## Print table
print(chl_bal_tab_opp)

### Table A3, Panel A ###
# Create row labels
arg_results_table <- data.frame(c("Gain frame", "Loss frame", "Difference in means", "$t$-statistic", "$p$-value"))

# Calculate t-test results
for(i in c("No opponent", "Chile", "UK")){
  # Gain frame observations
  gain_obs <- subset(arg_dat, frame_loss == 0 & opp_name == i)$gamble
  # Loss frame observations
  loss_obs <- subset(arg_dat, frame_loss == 1 & opp_name == i)$gamble
  # Gain frame mean
  gain_mean <- mean(gain_obs, na.rm = T)
  # Loss frame mean
  loss_mean <- mean(loss_obs, na.rm = T)
  # Subtract means
  diff_means <- loss_mean - gain_mean
  # Conduct t-test
  ttest_out <- t.test(loss_obs, gain_obs)
  # Bind results to table
  arg_results_table <- cbind(arg_results_table, round(c(gain_mean, loss_mean, diff_means, ttest_out$statistic, ttest_out$p.value), 3))
}

# Rename columns
colnames(arg_results_table) <- c("Quantity", "No opponent", "Chile", "UK")

# Print table
print(arg_results_table)


## Table A3, Panel B ##
# Create row labels
arg_results_within <- data.frame(c("Difference in means, gain frame", "$t$-statistic", "$p$-value", "Difference in means, loss frame", "$t$-statistic", "$p$-value"))

# Calculate t-test results: loop through every pairing of opponent conditions (Chile vs. none, UK vs. Chile, UK vs. none)
for(i in c(1:3)){
  # Gain frame responses, opponent condition 1
  obs1 <- subset(arg_dat, frame_loss == 0 & opp_name == combn(c("No opponent", "Chile", "UK"), 2)[1, i])$gamble
  # Gain frame responses, opponent condition 2
  obs2 <- subset(arg_dat, frame_loss == 0 & opp_name == combn(c("No opponent", "Chile", "UK"), 2)[2, i])$gamble
  # Difference in means, gain frame
  diff_means <- mean(obs2, na.rm = T) - mean(obs1, na.rm = T)
  # T-test, gain frame
  ttest_out <- t.test(obs2, obs1)
  # Loss frame responses, opponent condition 1
  obs1.loss <- subset(arg_dat, frame_loss == 1 & opp_name == combn(c("No opponent", "Chile", "UK"), 2)[1, i])$gamble
  # Loss frame responses, opponent condition 2
  obs2.loss <- subset(arg_dat, frame_loss == 1 & opp_name == combn(c("No opponent", "Chile", "UK"), 2)[2, i])$gamble
  # Difference in means, loss frame
  diff_means.loss <- mean(obs2.loss, na.rm = T) - mean(obs1.loss, na.rm = T)
  # T-test, loss frame
  ttest_out.loss <- t.test(obs2.loss, obs1.loss)
  # Bind to table
  arg_results_within <- cbind(arg_results_within, round(c(diff_means, ttest_out$statistic, ttest_out$p.value, diff_means.loss, ttest_out.loss$statistic, ttest_out.loss$p.value), 3))
}

# Rename columns
colnames(arg_results_within) <- c("Quantity", "Chile - None", "UK - None", "UK - Chile")

# Re-order columns
arg_results_within <- arg_results_within %>% dplyr::select(Quantity, "Chile - None", "UK - Chile", "UK - None")

# Print table
print(arg_results_within)


### Table A4, Panel A ###
# Create row labels
chl_results_table <- data.frame(c("\\% prefer risky option, gain frame", "\\% prefer risky option, loss frame", "Diff. in means (framing effect)", "$t$-statistic", "$p$-value"))

# Calculate t-test results
for(i in c("Argentina", "UK")){
  # Gain frame responses
  gain_obs <- subset(chl_dat, frame_loss == 0 & opp_name == i)$gamble
  # Loss frame responses
  loss_obs <- subset(chl_dat, frame_loss == 1 & opp_name == i)$gamble
  # Gain frame mean
  gain_mean <- mean(gain_obs, na.rm = T)
  # Loss frame mean
  loss_mean <- mean(loss_obs, na.rm = T)
  # Difference in means
  diff_means <- loss_mean - gain_mean
  # T-test
  ttest_out <- t.test(loss_obs, gain_obs)
  # Bind results to table
  chl_results_table <- cbind(chl_results_table, round(c(gain_mean, loss_mean, diff_means, ttest_out$statistic, ttest_out$p.value), 3))
}

# Rename columns
colnames(chl_results_table) <- c("Quantity", "Argentina", "UK")

# Print table
print(chl_results_table)


### Table A4, Panel B ###
chl_results_within <- data.frame(Quantity = c("Difference in means, gain frame", "$t$-statistic", "$p$-value", "Difference in means, loss frame", "$t$-statistic", "$p$-value"))

# Gain frame responses, Argentina opponent
obs1 <- subset(chl_dat, frame_loss == 0 & opp_name == "Argentina")$gamble
# Gain frame responses, UK opponent
obs2 <- subset(chl_dat, frame_loss == 0 & opp_name == "UK")$gamble
# Difference in means
diff_means <- mean(obs2, na.rm = T) - mean(obs1, na.rm = T)
# T-test
ttest_out <- t.test(obs2, obs1)
# Loss frame responses, Argentina opponent
obs1.loss <- subset(chl_dat, frame_loss == 1 & opp_name == "Argentina")$gamble
# Loss frame responses, UK opponent
obs2.loss <- subset(chl_dat, frame_loss == 1 & opp_name == "UK")$gamble
# Difference in means
diff_means.loss <- mean(obs2.loss, na.rm = T) - mean(obs1.loss, na.rm = T)
# T-test
ttest_out.loss <- t.test(obs2.loss, obs1.loss)

# Bind results to table
chl_results_within$`UK - Argentina` <- round(c(diff_means, ttest_out$statistic, ttest_out$p.value, diff_means.loss, ttest_out.loss$statistic, ttest_out.loss$p.value), 3)

# Print table
print(chl_results_within)


### Table A5 ###
list(lm(gamble ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk), data = arg_dat),
     lm(gamble ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat),
     glm(gamble ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat, family = binomial(link = "logit"))) %>%
stargazer(table.placement = "!t",
          type = "text",
          title="Argentina framing experiment main results", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:results_arg", 
          model.names = F,
          dep.var.labels = "Prefer risky policy option",
          column.labels = c("\\textit{OLS}", "\\textit{Logit}"), 
          column.separate = c(2, 1),
          covariate.labels=c("Loss frame", "Chile", "UK", "Loss frame x Chile","Loss frame x UK", "Age", "Female", "White", "Higher ed.", "Buenos Aires resid.", "Ideology","National attach.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table A6 ###
list(lm(gamble ~ frame_loss + opp_uk + I(frame_loss*opp_uk), data = chl_dat %>% filter(!is.na(opp_uk))),
     lm(gamble ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat %>% filter(!is.na(opp_uk))),
     glm(gamble ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat %>% filter(!is.na(opp_uk)), family = binomial(link = "logit"))) %>%
stargazer(table.placement = "!t",
          type = "text",
          title="Chile framing experiment main results", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:chl_results", 
          model.names = F,
          dep.var.labels = "Prefer risky policy option",
          column.labels = c("\\textit{OLS}", "\\textit{Logit}"), 
          column.separate = c(2, 1),
          covariate.labels=c("Loss frame", "UK", "Loss frame x UK", "Age", "Female", "White", "Higher educ.", "Santiago", "Ideology", "National identif.", "Militant assert.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table B1 ###
## Create table
chl_bal_tab_hist <- chl_dat %>%
  # Group by treatment cells
  group_by(frame_loss, historical, hi_val) %>%
  # Get mean value of each demographic variable by treatment cell
  summarise(across(c(age, female, race_white, educ_superior, santiago, ideology_scale), ~ mean(.x, na.rm = T))) %>%
  ungroup() %>% dplyr::select(age:ideology_scale) %>% pivot_longer(cols = age:ideology_scale) %>% group_by(name) %>%
  # Get minimum and maximum mean values of each demographic variable
  summarise(min = min(value), max = max(value)) %>%
  # Re-order table
  mutate(name = factor(name, levels = c("age", "female", "race_white", "educ_superior", "santiago", "ideology_scale"))) %>% arrange(name) %>%
  # Merge population means of each demographic variable
  mutate(target_means = c(42.38, # Age; LAPOP 2021 
                          8946991/17520239, # Female; Census 2017
                          0.437, # White; LAPOP 2021
                          0.298, # Superior educ.; Census 2017
                          7112808/17520239, # Santiago residents as share of population; Census 2017
                          5.054), # Ideology; LAPOP 2018
         # Merge overall sample mean of each demographic variable
         sample_means = chl_dat %>% filter(!is.na(historical)) %>% dplyr::select(age, female, race_white, educ_superior, santiago, ideology_scale) %>% colMeans(na.rm = T))

## Calculate standardized mean differences using 'cobalt' package
# Drop all observations with missing demographic data to use in cobalt
chl_dat_sub_hist <- chl_dat %>% dplyr::select(frame_loss, historical, hi_val, gamble, age, female, race_white, educ_superior, santiago, ideology_scale) %>% na.omit

# Estimate standardized mean differences using cobalt
chl_cobalt_hist <- bal.tab(chl_dat_sub_hist %>% dplyr::select(age:ideology_scale), treat = paste0(chl_dat_sub_hist$frame_loss, chl_dat_sub_hist$historical, chl_dat_sub_hist$hi_val))

# Incorporate standardized mean difference estimates from cobalt into table
chl_bal_tab_hist$std_mean_diff <- chl_cobalt_hist$Balance.Across.Pairs[,2] %>% round(3)

## Clean up 
# Reorder columns
chl_bal_tab_hist <- dplyr::select(chl_bal_tab_hist, name, target_means, sample_means, max, min, std_mean_diff)

# Label rows and columns
chl_bal_tab_hist$name <- c("Age", "Female", "White", "Superior ed.", "Santiago", "Right-wing ideology (1-10)")
colnames(chl_bal_tab_hist) <- c("Covariate", "Pop. mean", "Sample mean", "Max. cell mean", "Min. cell mean", "Std. mean diff.")

## Print table
print(chl_bal_tab_hist)


### Table B2 ###
list(lm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val), data = chl_dat), 
     lm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data= chl_dat), 
     glm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "logit"), data = chl_dat)) %>%
stargazer(type = "text",
          table.placement = "H",
          title="Historical framing experiment main results", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:results_hist_chile", 
          model.names = F,
          dep.var.labels = c("Prefer risky policy option"),
          column.labels = c("\\textit{OLS}", "\\textit{Logit}"), 
          column.separate = c(2, 1),
          covariate.labels=c("Loss frame", "Historical claim", "High value", 
                             "Loss frame x Hist. claim", "Loss frame x High value", "Hist. claim x High value",
                             "Loss frame x Hist. claim x High value",
                             "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", 
                             "National identif.", "Militant assert.",
                             "International trust", "Social trust", 
                             "Intercept"),
          no.space = T)

# F-test for historical ownership and economic value treatments, Model 1
anova(lm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val), data = chl_dat), 
      lm(gamble ~ frame_loss, data = chl_dat %>% filter(!is.na(historical), !is.na(hi_val))))

# F-test for historical ownership and economic value treatments, Model 2
anova(lm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat), 
      lm(gamble ~ frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat %>% filter(!is.na(historical), !is.na(hi_val))))

# Likelihood ratio test for historical ownership and economic value treatments, Model 3
anova(glm(gamble ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "logit"), data = chl_dat),
      glm(gamble ~ frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "logit"), data = chl_dat %>% filter(!is.na(historical), !is.na(hi_val))),
      test = "Chisq")


### Table B3 ###
list(lm(perceive_own_bin ~ historical + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     lm(perceive_own_bin ~ historical + hi_val + I(historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     lm(perceive_own_bin ~ historical + hi_val + frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     glm(perceive_own_bin ~ historical + hi_val + I(historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data= chl_dat, family = binomial(link = "logit"))) %>%
stargazer(type = "text",
          title="Effect of treatments on perceived ownership", 
          label="tab:treat_ownership", 
          table.placement = "!t",
          omit.stat = c("f", "ser", "adj.rsq", "ll"),
          dep.var.labels = c("Perceived ownership"),
          model.names = F,
          column.labels = c("\\textit{OLS}", "\\textit{Logit}"),
          column.separate = c(3, 1),
          covariate.labels=c("Historical frame", "High value", "Historical x High value", "Loss frame", "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", "National identif.", "Militant assert.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table B4 ###
## IV regression, OLS (Model 3)
hist_itt <- ivreg(gamble ~ perceive_own_bin + frame_loss + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
                    historical + frame_loss + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat)

## IV probit regression (Model 4)
# First stage
ivprobit_s1 <- lm(perceive_own_bin ~ historical + frame_loss + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat)

# Generate predicted values of perceived ownership: Copy data
chl_dat_probit <- chl_dat
# Replace observed values of perceived ownership
chl_dat_probit$perceive_own_bin <- NA 
# Insert predicted values of perceived ownership
chl_dat_probit$perceive_own_bin[as.numeric(labels(predict(ivprobit_s1)))] <- predict(ivprobit_s1) 

# Second stage
ivprobit_s2 <- glm(gamble ~ perceive_own_bin + frame_loss + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, family = binomial(link = "probit"), data = chl_dat_probit)

## Create table
list(lm(gamble ~ perceive_own_bin, data = subset(chl_dat, !is.na(historical))), 
     lm(gamble ~ perceive_own_bin + historical + frame_loss + hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = subset(chl_dat, !is.na(historical))), 
     hist_itt,
     ivprobit_s2) %>% 
stargazer(type = "text",
          table.placement = "H",
          title="Effect of perceived ownership on risk preferences", 
          omit.stat=c("LL", "ser", "f", "adj.rsq", "rsq", "aic"),
          digits=3, 
          label="tab:ownership_risk", 
          model.names = F,
          dep.var.labels = c("Prefer risky policy option"),
          column.labels = c("\\textit{OLS}", "\\textit{IV}", "\\textit{IV Probit}"), 
          column.separate = c(2, 1, 1),
          covariate.labels=c("Perceived ownership", "Historical ownership", "Loss frame", "High value", "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", "National identif.", "Militant assert.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table C1 ###
## Panel A
# Estimate models
arg_mil <- glm(policy_mil ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat, family = binomial)
arg_ic <- glm(policy_ic ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat, family = binomial)
arg_reneg <- glm(policy_reneg ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat, family = binomial)
arg_wait <- glm(policy_wait ~ frame_loss + chile + uk + I(frame_loss*chile) + I(frame_loss*uk) + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust, data = arg_dat, family = binomial)

# Create table of likelihood ratio tests
arg_ftests <- matrix(nrow = 3, ncol = 4)

# Populate table by looping through each model and treatment variable
for(i in 1:4){
  # Pull model
  full_mod <- list(arg_mil, arg_ic, arg_reneg, arg_wait)[[i]]
  
  # Extract dependent variable
  dv <- as.character(full_mod$formula)[2]
  
  for(j in 1:3){
    # Pull treatment variable to exclude
    exclude_treat <- c("frame_loss", "chile", "uk")[j]
    
    # Remove any term that includes the treatment variable
    include_treats <- c("frame_loss", "chile", "uk", "I(frame_loss*chile)", "I(frame_loss*uk)")[!grepl(exclude_treat, c("frame_loss", "chile", "uk", "I(frame_loss*chile)", "I(frame_loss*uk)"))]
    
    # Estimate parsimonious model with controls
    compare_mod <- glm(as.formula(paste0(dv, "~", paste(c(include_treats, "age", "female", "race_white", "university", "buenos_aires", "ideology", "natattach", "int_trust", "soc_trust"), collapse = "+"))), data = arg_dat, family = binomial)
   
    # Conduct likelihood ratio test between the full model and parsimonious model
    ftest <- anova(full_mod, compare_mod, test = "Chisq")
    
    # Extract p-value
    p <- formatC(ftest$`Pr(>Chi)`[2], digits = 3, format = "f")
    
    # Attach significance stars and insert into table
    if(p < 0.1 & p >= 0.05){
      arg_ftests[j, i] <- paste0(p, "*")
    }else if(p < 0.05 & p >= 0.01){
      arg_ftests[j, i] <- paste0(p, "**")
    }else if(p < 0.01){
      arg_ftests[j, i] <- paste0(p, "***")
    }else {
      arg_ftests[j, i] <- p
    }
  }
}

# Label rows and columns
rownames(arg_ftests) <- c("Loss frame", "Chile", "UK")
colnames(arg_ftests) <- c("Military action", "Int'l court", "Renegotiate", "Wait")

# Print table
print(arg_ftests)


## Panel B
# Estimate models
opp_mil <- glm(policy_mil ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
opp_ic <- glm(policy_ic ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
opp_reneg <- glm(policy_reneg ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
opp_wait <- glm(policy_wait ~ frame_loss + opp_uk + I(frame_loss*opp_uk) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)

# Create table of likelihood ratio tests
opp_ftests <- matrix(nrow = 2, ncol = 4)

# Populate table by looping through each model and treatment variable
for(i in 1:4){
  # Pull model
  full_mod <- list(opp_mil, opp_ic, opp_reneg, opp_wait)[[i]]
  
  # Extract dependent variable
  dv <- as.character(full_mod$formula)[2]
  
  for(j in 1:2){
    # Pull treatment variable to exclude
    exclude_treat <- c("frame_loss", "opp_uk")[j]
    
    # Remove any term that includes the treatment variable
    include_treats <- c("frame_loss", "opp_uk", "I(frame_loss*opp_uk)")[!grepl(exclude_treat, c("frame_loss", "opp_uk", "I(frame_loss*opp_uk)"))]
    
    # Estimate parsimonious model with controls
    compare_mod <- glm(as.formula(paste0(dv, "~", paste(c(include_treats, "age", "female", "race_white", "educ_superior", "santiago", "ideology", "natid_scale", "mi_scale", "int_trust", "soc_trust"), collapse = "+"))), data = subset(chl_dat, !is.na(opp_uk)), family = binomial)
    
    # Conduct likelihood ratio test between the full model and parsimonious model
    ftest <- anova(full_mod, compare_mod, test = "Chisq")
    
    # Extract p-value
    p <- formatC(ftest$`Pr(>Chi)`[2], digits = 3, format = "f")
    
    # Attach significance stars and insert into table
    if(p < 0.1 & p >= 0.05){
      opp_ftests[j, i] <- paste0(p, "*")
    }else if(p < 0.05 & p >= 0.01){
      opp_ftests[j, i] <- paste0(p, "**")
    }else if(p < 0.01){
      opp_ftests[j, i] <- paste0(p, "***")
    }else {
      opp_ftests[j, i] <- p
    }
  }
}

# Label rows and columns
rownames(opp_ftests) <- c("Loss frame", "UK")
colnames(opp_ftests) <- c("Military action", "Int'l court", "Renegotiate", "Wait")

# Print table
print(opp_ftests)


## Panel C
# Estimate models
hist_mil <- glm(policy_mil ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
hist_ic <- glm(policy_ic ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
hist_reneg <- glm(policy_reneg ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)
hist_wait <- glm(policy_wait ~ frame_loss + historical + hi_val + I(frame_loss*historical) + I(frame_loss*hi_val) + I(historical*hi_val) + I(frame_loss*historical*hi_val) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial)

# Create table of likelihood ratio tests
hist_ftests <- matrix(nrow = 3, ncol = 4)

# Populate table by looping through each model and treatment variable
for(i in 1:4){
  # Pull model
  full_mod <- list(hist_mil, hist_ic, hist_reneg, hist_wait)[[i]]
  
  # Extract dependent variable
  dv <- as.character(full_mod$formula)[2]
  
  for(j in 1:3){
    # Pull treatment variable to exclude
    exclude_treat <- c("frame_loss", "historical", "hi_val")[j]
    
    # Remove any term that includes the treatment variable
    include_treats <- c("frame_loss", "historical", "hi_val", "I(frame_loss*historical)", "I(frame_loss*hi_val)", "I(historical*hi_val)", "I(frame_loss*historical*hi_val)")[!grepl(exclude_treat, c("frame_loss", "historical", "hi_val", "I(frame_loss*historical)", "I(frame_loss*hi_val)", "I(historical*hi_val)", "I(frame_loss*historical*hi_val)"))]
    
    # Estimate parsimonious model with controls
    compare_mod <- glm(as.formula(paste0(dv, "~", paste(c(include_treats, "age", "female", "race_white", "educ_superior", "santiago", "ideology", "natid_scale", "mi_scale", "int_trust", "soc_trust"), collapse = "+"))), data = subset(chl_dat, !is.na(historical)), family = binomial)
    
    # Conduct likelihood ratio test between the full model and parsimonious model
    ftest <- anova(full_mod, compare_mod, test = "Chisq")
    
    # Extract p-value
    p <- formatC(ftest$`Pr(>Chi)`[2], digits = 3, format = "f")
    
    # Attach significance stars and insert into table
    if(p < 0.1 & p >= 0.05){
      hist_ftests[j, i] <- paste0(p, "*")
    }else if(p < 0.05 & p >= 0.01){
      hist_ftests[j, i] <- paste0(p, "**")
    }else if(p < 0.01){
      hist_ftests[j, i] <- paste0(p, "***")
    }else {
      hist_ftests[j, i] <- p
    }
  }
}

# Label rows and columns
rownames(hist_ftests) <- c("Loss frame", "Historical", "High value")
colnames(hist_ftests) <- c("Military action", "Int'l court", "Renegotiate", "Wait")

# Print table
print(hist_ftests)


### Table C2 ###
stargazer(arg_mil, arg_ic, arg_reneg, arg_wait,
          type = "text",
          table.placement = "!t",
          title="Treatment effects on perception of risky policy option, Argentina (logit)", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:gamble_policy_arg", 
          dep.var.labels = c("Military action", "Int'l court", "Renegotiate", "Wait"),
          covariate.labels=c("Loss frame", "Chile", "UK", "Loss frame x Chile","Loss frame x UK", "Age", "Female", "White", "Higher ed.", "Buenos Aires", "Ideology", "National attachment", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table C3 ###
stargazer(opp_mil, opp_ic, opp_reneg, opp_wait,
          type = "text",
          table.placement = "!t",
          title="Experiment 1 treatment effects on perception of risky policy option, Chile (logit)", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:gamble_policy_chile1", 
          dep.var.labels = c("Military action", "Int'l court", "Renegotiate", "Wait"),
          covariate.labels=c("Loss frame", "UK", "Loss frame x UK", "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", "National identification", "Militant assert.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table C4 ###
stargazer(hist_mil, hist_ic, hist_reneg, hist_wait,
          type = "text",
          table.placement = "!t",
          title="Experiment 2 treatment effects on perception of risky policy option, Chile (logit)", 
          omit.stat=c("LL", "ser", "f", "adj.rsq"),
          digits=3, 
          label="tab:gamble_policy_chile2", 
          dep.var.labels = c("Military action", "Int'l court", "Renegotiate", "Wait"),
          covariate.labels=c("Loss frame", "Historical", "High value", "Loss frame x Historical", "Loss frame x High value", "Historical x High value", "Loss frame x Historical x High value", "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology", "National identification", "Militant assert.", "International trust", "Social trust", "Intercept"),
          no.space = T)


### Table C5 ###
# Function to examine whether framing effect is mediated by policy interpretation
frame_mediate <- function(country = "argentina", opponent = "uk", policy = "mil", sim = 1000){
  
  # Set up data and regression inputs depending on country
  if(country == "argentina"){
    dat <- arg_dat
    covars <- c("age", "female", "race_white", "university", "buenos_aires", "ideology", "natattach", "int_trust", "soc_trust")
    treats <- "frame_loss*chile + frame_loss*uk"
  }else if(country == "chile"){
    dat <- subset(chl_dat, !is.na(opp_uk))
    covars <- c("age", "female", "race_white", "educ_superior", "santiago", "ideology", "natid_scale", "mi_scale", "int_trust", "soc_trust")
    treats <- "frame_loss*opp_uk"
  }
  
  # Policy interpretation variable
  policy_var <- paste0("policy_", policy)
  
  # Mediating model
  med.mod <- glm(as.formula(paste(c(paste(policy_var, "~", treats), covars), collapse = "+")), data = subset(dat, !is.na(gamble)), family = "binomial")
  
  # Outcome model
  out.mod <-  glm(as.formula(paste(c(paste("gamble ~", policy_var, "+", treats), covars), collapse = "+")), data = dat, family = "binomial")
  
  # Conduct mediation analysis depending on country and opponent condition
  if(country == "argentina"){
    if(opponent == "none"){
      med.out.frame.effect <- mediate(med.mod, out.mod, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(chile = 0, uk = 0))
    }else if(opponent == "chile"){
      med.out.frame.effect <- mediate(med.mod, out.mod, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(chile = 1, uk = 0))
    }else if(opponent == "uk"){
      med.out.frame.effect <- mediate(med.mod, out.mod, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(chile = 0, uk = 1))
    }
  } else if(country == "chile"){
    if(opponent == "argentina"){
      med.out.frame.effect <- mediate(med.mod, out.mod, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(opp_uk = 0))
    }else if(opponent == "uk"){
      med.out.frame.effect <- mediate(med.mod, out.mod, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(opp_uk = 1))
    }
  }
  return(med.out.frame.effect)
}

# Function to extract p-values and proportion mediated from mediation analysis
frame_med_prop <- function(pol, ctry, opp, sims = 1000){
  # Mediation analysis
  med.obj <- frame_mediate(country = ctry, opponent = opp, policy = pol, sim = sims)
  
  # P-value
  p <- formatC(med.obj$n.avg.p, digits = 3, format = "f")
  
  # Proportion mediated
  prop <- formatC(med.obj$n.avg, digits = 3, format = "f")
  
  # Add stars of significance
  if(p < 0.1 & p >= 0.05){
    return(c(paste0(prop, "*"), paste0("(",p, ")")))
  }else if(p < 0.05 & p>= 0.01){
    return(c(paste0(prop, "**"), paste0("(",p, ")")))
  }else if(p < 0.01){
    return(c(paste0(prop, "***"), paste0("(",p, ")")))
  }else{
    return(c(prop, paste0("(",p, ")")))
  }
}

# Function to cycle through policy interpretations and conduct mediation analysis for each
frame_med_prop_opp <- function(opponent, country, sim){
  return(sapply(c("mil", "ic", "reneg", "wait"), frame_med_prop, ctry = country, opp = opponent, sims = sim))
}

# Function to conduct mediation analysis for historical ownership experiment
frame_mediate_hist <- function(policy, sim = 1000, hist = NULL, val = NULL){
  # Policy interpretation variable
  policy_var <- paste0("policy_", policy)
  
  # Mediating model
  med.hist <- glm(as.formula(paste(policy_var, "~ frame_loss*historical*hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust")), data = subset(chl_dat, !is.na(gamble) & !is.na(historical)), family = "binomial")
  
  # Outcome model
  out.hist <- glm(as.formula(paste("gamble ~", policy_var, "+ frame_loss*historical*hi_val + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust")), data = subset(chl_dat, !is.na(historical)), family = "binomial")
  
  # Conduct mediation analysis
  if(is.null(hist) & is.null(val)){
    med.obj <- mediate(med.hist, out.hist, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim)
  } else {
    med.obj <- mediate(med.hist, out.hist, treat = "frame_loss", mediator = policy_var, robustSE = T, sims = sim, covariates = list(historical = hist, hi_val = val))                 
  }
  
  # P-value
  p <- formatC(med.obj$n.avg.p, digits = 3, format = "f")
  
  # Proportion mediated
  prop <- formatC(med.obj$n.avg, digits = 3, format = "f")
  
  # Add stars of significance
  if(p < 0.1 & p >= 0.05){
    return(c(paste0(prop, "*"), paste0("(",p, ")")))
  }else if(p < 0.05 & p>= 0.01){
    return(c(paste0(prop, "**"), paste0("(",p, ")")))
  }else if(p < 0.01){
    return(c(paste0(prop, "***"), paste0("(",p, ")")))
  }else{
    return(c(prop, paste0("(",p, ")")))
  }
}

# Set seed
set.seed(1234)

# Construct table
frame_med <- cbind(# Argentina sample
  sapply(c("none", "chile", "uk"), frame_med_prop_opp, "argentina", 1000),
  # Chile sample
  sapply(c("argentina", "uk"), frame_med_prop_opp, "chile", 1000),
  # Historical ownership experiment
  as.vector(sapply(c("mil", "ic", "reneg", "wait"), frame_mediate_hist, sim = 1000)))

# Label columns and rows
colnames(frame_med) <- c("No opponent", "Chile", "UK", "Argentina", "UK", "Pooled")
rownames(frame_med) <- c("Military action", "", "International court", "", "Renegotiate", "", "Wait", "")

# Print table
print(frame_med)


### Table C6 ###
# Function to examine whether policy interpretation mediates opponent effects
opponent_mediate <- function(policy = "mil", country = "argentina", opponent = "uk", frame = "gain", sim = 1000){
  # Policy interpretation variable
  policy_var <- paste0("policy_", policy)
  
  # Mediating and outcome models
  if(country == "argentina"){
    med.mod <- glm(as.formula(paste(policy_var, "~ frame_loss*chile + frame_loss*uk + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust")), data = subset(arg_dat, !is.na(gamble)), family = "binomial")
    out.mod <- glm(as.formula(paste("gamble ~", policy_var, "+", "frame_loss*chile + frame_loss*uk + age + female + race_white + university + buenos_aires + ideology + natattach + int_trust + soc_trust")), data = arg_dat, family = "binomial")
    med.out.opp.effect <- mediate(med.mod, out.mod, treat = opponent, mediator = policy_var, robustSE = T, sims = sim, covariates = list(frame_loss = ifelse(frame == "loss", 1, 0)))
  }else{
    med.mod <- glm(as.formula(paste(policy_var, "~ frame_loss*opp_uk + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust")), data = subset(chl_dat, !is.na(gamble)), family = "binomial")
    out.mod <- glm(as.formula(paste("gamble ~", policy_var, "+", "frame_loss*opp_uk + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust")), data = chl_dat, family = "binomial")
    med.out.opp.effect <- mediate(med.mod, out.mod, treat = "opp_uk", mediator = policy_var, robustSE = T, sims = sim, covariates = list(frame_loss = ifelse(frame == "loss", 1, 0)))
  }
  return(med.out.opp.effect)
}

# Function to extract p-value and proportion mediated
opp_med_prop <- function(pol, ctry, opp, f, sims = 1000){
  # Create mediation object
  med.obj <- opponent_mediate(country = ctry, opponent = opp, policy = pol, frame = f, sim = sims)
  
  # Extract p-value
  p <-  formatC(med.obj$n.avg.p, digits = 3, format = "f")
  
  # Extract proportion mediated
  prop <- formatC(med.obj$n.avg, digits = 3, format = "f")
  
  # Attach stars of significance
  if(p < 0.1 & p >= 0.05){
    return(c(paste0(prop, "*"), paste0("(",p, ")")))
  }else if(p < 0.05 & p>= 0.01){
    return(c(paste0(prop, "**"), paste0("(",p, ")")))
  }else if(p < 0.01){
    return(c(paste0(prop, "***"), paste0("(",p, ")")))
  }else{
    return(c(prop, paste0("(",p, ")")))
  }
}

# Set seed
set.seed(1234)

# Construct table
opp_med <- cbind(sapply(c("chile", "uk"), function(opponent){sapply(c("mil", "ic", "reneg", "wait"), opp_med_prop, ctry = "argentina", opp = opponent, f = "gain", sims = 1000)}),
                 sapply(c("gain", "loss"), function(frame){sapply(c("mil", "ic", "reneg", "wait"), opp_med_prop, ctry = "chile", f = frame, sims = 1000)}))

# Label columns and rows
colnames(opp_med) <- c("Chile: gain", "UK: gain", "UK: gain", "UK: loss")
rownames(opp_med) <- c("Military action", "", "International court", "", "Renegotiate", "", "Wait", "")

# Print table
print(opp_med)


### Table D1 ###
# Calculate propensity to interpret risky option as military action using pre-treatment covariates
mil_prop_mod <- glm(policy_mil ~ frame_loss + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat, family = binomial(link = "logit"))
mil_prop <- predict(mil_prop_mod, type = "response")
chl_dat$mil_prop <- NA
chl_dat$mil_prop[as.integer(labels(mil_prop))] <- mil_prop

# Create table
list(ivreg(force_ab ~ gamble + policy_mil + I(gamble * policy_mil) + I(gamble * mi_scale) + I(gamble * natid_scale) + I(gamble * age) + I(gamble * female) + I(gamble * int_trust) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
             frame_loss + policy_mil + I(frame_loss * policy_mil) + I(frame_loss * mi_scale) + I(frame_loss * natid_scale) + I(frame_loss * age) + I(frame_loss * female) + I(frame_loss * int_trust) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat),
     ivreg(force_ab ~ gamble + policy_mil + I(gamble * policy_mil) + mil_prop + I(gamble * mil_prop) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust | 
             frame_loss + policy_mil + I(frame_loss * policy_mil) + mil_prop + I(frame_loss * mil_prop) + age + female + race_white + educ_superior + santiago + ideology + natid_scale + mi_scale + int_trust + soc_trust, data = chl_dat)) %>%
stargazer(type = "text",
          title = "Risk acceptance on use of force, supplemental results",
          label = "tab:force_supp",
          omit.stat = c("adj.rsq", "ser", "f"),
          dep.var.labels = "Support use of force",
          column.labels = c(rep("\\textit{IV}", 2)),
          model.names = F,
          order = c("^gamble", "frame_loss", "^policy", "gamble"),
          covariate.labels=c("Risky option", "Mil. action",
                             "Risky option x Mil. action", "Risky option x Mil. assert.", 
                             "Risky option x Nat'l iden.", "Risky option x Age", "Risky option x Female", "Risky option x Int'l trust", "Risky option x propensity mil. action",
                             "Propensity mil. action",
                             "Age", "Female", "White", "Higher ed.", "Santiago", "Ideology",
                             "National identif.", "Mil. assertiveness",
                             "Int'l trust", "Social trust"),
          no.space = T)

### Figure E1 ###
# Create list to store marginal effect outputs
emo_list <- list()

# Estimate marginal opponent effect as function of how differently the respondent feels toward the UK vs. Argentina
for(frame in c(0,1)){ # Loop through each frame
  for(emo in 1:4){# Loop through each emotion
    df <- interflex(estimator = "linear", method = "logit", 
                    # Subset sample to one frame at a time
                    subset(chl_dat, frame_loss == frame), 
                    
                    # Outcome: risk acceptance
                    Y = "gamble", 
                    
                    # Treatment: whether the opponent was the UK (opp_all = 1) or Argentina (opp_all = 0)
                    D = "opp_all", 
                    
                    # Moderator: difference between emotion rating for UK vs. Argentina
                    X = c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[emo],
                    
                    # Controls: each of the other emotions, their interactions with the opponent treatment, and demographic and dispositional covariates
                    Z = c(c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[-emo], paste0(c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[-emo], "xopp_all"), "age", "female", "race_white", "educ_superior", "santiago", "ideology", "natid_scale", "mi_scale", "int_trust", "soc_trust"), 
                    na.rm = T, treat.type = "discrete", base = 0, figure = F)$est.lin[[1]] %>% data.frame %>%
     
      # Label emotion and frame in output
       mutate(emotion = c("Warmth", "Fear", "Anger", "Confidence")[emo], 
              frame = c("Gain", "Loss")[frame + 1])
    
    # Store output into list
    emo_list <- c(emo_list, list(df))
  }
}

# Bind marginal effect output list into a single table and reorder emotions for plotting
emo_df <- do.call(rbind, emo_list) %>% mutate(emotion = factor(emotion, levels = unique(emotion)))

# Create plot
ggplot(emo_df, aes(x = X, y = ME, ymin = lower.CI.95.., ymax = upper.CI.95..)) + 
  geom_line() + geom_ribbon(alpha = 0.3, color = NA) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  xlab("Relative emotion toward UK (UK - Argentina)") + ylab("Marginal effect of UK as opponent") +
  facet_grid(frame~emotion, scales = "free_x")


### Figure E2 ###
# Create interaction between framing and opponent treatments
chl_dat$frame_lossxopp_all <- chl_dat$frame_loss*chl_dat$opp_all

# For each emotion, create two-way interaction with framing and three-way interaction with framing and opponent treatments
for(e in c("warm_diff", "fear_diff", "anger_diff", "strong_diff")){
  chl_dat[[paste0(e, "xframe_loss")]] <- chl_dat[[e]]*chl_dat[["frame_loss"]]
  chl_dat[[paste0(e, "xframe_lossxopp_all")]] <- chl_dat[[e]]*chl_dat[["frame_loss"]]*chl_dat[["opp_all"]]
}

# Create list to store marginal effect outputs
emo_frame <- list()

# Estimate marginal framing effect as function of how differently the respondent feels toward the UK vs. Argentina
for(emo in 1:4){ # Loop through each emotion
  df <- interflex(estimator = "linear", method = "logit", chl_dat, 
                  # Outcome: Risk acceptance
                  Y = "gamble", 
                 
                  # Treatment: Difference in framing effect depending on whether the opponent is the UK vs. Argentina
                  D = "frame_lossxopp_all", 
                  
                  # Moderator: difference in emotion rating for UK vs. Argentina
                  X = c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[emo],
                  
                  # Controls: each uninteracted treatment variable, each of the other emotions, their two-way interactions with the framing treatment and opponent treatment, their three-way interactions with the framing and opponent treatments, and demographic and dispositional covariates
                  Z = c("frame_loss", "opp_all", c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[-emo], paste0(c("warm_diff", "fear_diff", "anger_diff", "strong_diff"), "xframe_loss"), paste0(c("warm_diff", "fear_diff", "anger_diff", "strong_diff"), "xopp_all"), paste0(c("warm_diff", "fear_diff", "anger_diff", "strong_diff")[-emo], "xframe_lossxopp_all"), "age", "female", "race_white", "educ_superior", "santiago", "ideology", "natid_scale", "mi_scale", "int_trust", "soc_trust"), 
                  na.rm = T, treat.type = "discrete", base = 0, figure = F)$est.lin[[1]] %>% data.frame %>%
    
    # Label emotion in output
    mutate(emotion = c("Warmth", "Fear", "Anger", "Confidence")[emo])
  
  # Store output into list
  emo_frame <- c(emo_frame, list(df))
}

# Bind marginal effect output list into a single table and reorder emotions for plotting
emo_frame_df <- do.call(rbind, emo_frame) %>% mutate(emotion = factor(emotion, levels = unique(emotion)))

# Create plot
ggplot(emo_frame_df, aes(x = X, y = ME, ymin = lower.CI.95.., ymax = upper.CI.95..)) + 
  geom_line() + geom_ribbon(alpha = 0.3, color = NA) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  xlab("Relative emotion toward UK (UK - Argentina)") + ylab("Change in framing effect") +
  facet_grid(~emotion, scales = "free")


### Figure F1 ###
# Load data: Dispute data from Frederick et al. 2017, https://www.paulhensel.org/Data/ICOWprov.zip (provisional version 1.01); contiguity data from Hensel 2017, https://correlatesofwar.org/data-sets/direct-contiguity/ (v. 3.20)
icow_plot <- read.csv("icow_plot.csv")

## Left-hand panel
ggplot(icow_plot) + 
  geom_line(aes(x = year, y = begin_contig, color = "Started"), alpha = 0.4) +
  geom_line(aes(x = year, y = begin_smooth, color = "Started"), linewidth = 1.5) +
  geom_line(aes(x = year, y = end_contig, color = "Ended"), alpha = 0.4) +
  geom_line(aes(x = year, y = end_smooth, color = "Ended"), linewidth = 1.5)+
  geom_vline(aes(xintercept = 1945), linetype = "dotted")+
  xlab("Year") + ylab("Disputes per contiguous dyad")+
  theme(legend.title=element_blank(), plot.title=element_text(hjust=0.5), legend.position = "bottom")

## Right-hand panel
ggplot(icow_plot)+
  geom_line(aes(x = year, y = med.dur, color = "Median duration, ongoing claims"))+
  geom_line(aes(x = year, y = avg.end.dur, color = "Average duration, claims ended in preceding decade"))+
  geom_vline(aes(xintercept = 1945), linetype = "dotted")+
  xlab("Year")+
  ylab("Dispute duration")+
  theme(legend.title=element_blank(), plot.title=element_text(hjust=0.5), legend.position = "bottom", legend.direction = "vertical")


### Figure G1 ###
# Create table listing each moderating variable and its corresponding label
moderators_name <- data.frame(var = c("soc_trust", "int_trust", "authoritarian", "sdo", "ideology", "female", "natattach", "university"), 
                              moderator = c("Social trust", "International trust", "Authoritarianism", "SDO", "Ideology", "Female", "National attach.", "University"))

# Create function to generate predicted probabilities as function of moderator
het_arg_frame <- function(moderator){
  
  # Estimate logit regression
  reg_moderators <- glm(as.formula(paste(
    # Outcome: risk acceptance
    "gamble ~ ", 
    # Full interactions between moderator, framing, and opponent treatments
    paste(paste0(moderator, "*", c("frame_loss*chile", "frame_loss*uk")), collapse = "+"),
    # Demographic controls excluding ideology to avoid post-treatment bias
    "+ age + female + race_white + university + buenos_aires")), family = "binomial", data = arg_dat)
  
  ## Estimate predicted probabilities by framing treatment as function of moderator, separately for each opponent condition
  # Estimate at 0 and 1 if moderator is binary; estimate at 10 points if moderator is continuous
  if(moderator %in% c("female", "university")){
    x <- 2
  }else{
    x <- 10
  }
  margins <- c()
  for(opp in c("No opponent", "Chile", "UK")){ # Loop through opponent conditions
      for(f in 0:1){ # Loop through framing conditions
          # Estimate predicted probabilities
          marg <- data.frame(cplot(reg_moderators, x = moderator, what = "prediction", data = subset(arg_dat, opp_name == opp & frame_loss == f), draw = F, n = x))
          # Label opponent and frame
          marg$opponent <- opp
          marg$frame <- c("Gain", "Loss")[f+1]
          # Bind to results table
          margins <- rbind(margins, marg)
        }
      }
  
  # Add moderator name to output
  margins$var <- moderator
  return(margins)
}

# Apply function to every moderator
arg_marg_frame <- lapply(as.character(moderators_name$var), het_arg_frame) 
# Bind all results into a single table
arg_marg_frame <- bind_rows(arg_marg_frame, .id = "column_label")
# Merge results with moderator labels
arg_marg_frame <- merge(arg_marg_frame, moderators_name, by = "var", all.x = T)
# Re-order opponent conditions
arg_marg_frame$opponent <- factor(arg_marg_frame$opponent , levels = unique(arg_marg_frame$opponent))

# Create plot
ggplot()+
  geom_line(data = arg_marg_frame, aes(x = xvals, y = yvals, linetype = frame))+
  geom_ribbon(data = subset(arg_marg_frame, !var %in% c("female", "university")), aes(x = xvals, y = yvals, ymin = lower, ymax = upper, linetype = frame), alpha = 0.3) +
  geom_point(data = subset(arg_marg_frame, var %in% c("female", "university")), aes(x = xvals, y = yvals), size = 0.2)+
  geom_errorbar(data = subset(arg_marg_frame, var %in% c("female", "university")), aes(x = xvals, ymin = lower, ymax = upper), width = 0.07)+
  xlab("Moderator")+
  ylab("Predicted prob. of choosing risky option")+
  labs(fill = "", color = "", linetype = "")+
  theme(legend.position = "bottom")+
  facet_grid(moderator~opponent)


### Figure G2 ###
# Create table listing each moderating variable and its corresponding label
moderators_name_chl <- data.frame(var = c("soc_trust", "int_trust", "authoritarian", "ideology", "female", "natid_scale", "educ_superior"), 
                              moderator = c("Social trust", "International trust", "Authoritarianism", "Ideology", "Female", "National ident.", "Higher ed."))

# Create function to generate predicted probabilities as function of moderator
het_chl_opp_frame <- function(moderator){
  
  # Estimate logit regression
  reg_moderators <- glm(as.formula(paste(
    # Outcome: risk acceptance
    "gamble ~ ", 
    # Moderator fully interacted with treatment variables
    paste0(moderator, "*frame_loss*opp_uk"), 
    # Demographic controls excluding ideology to avoid post-treatment bias
    "+ age + female + race_white + educ_superior + santiago")), family = "binomial", data = chl_dat)
  
  ## Estimate predicted probabilities by framing treatment as function of moderator, separately for each opponent condition
  # Estimate at 0 and 1 if moderator is binary; estimate at 10 points if moderator is continuous
  if(moderator %in% c("female", "educ_superior")){
    x <- 2
  }else{
    x <- 10
  }
  margins <- c()
  for(opp in c("Argentina", "UK")){ # Loop through opponent conditions
    for(f in 0:1){ # Loop through framing conditions
      # Estimate predicted probabilities
      marg <- data.frame(cplot(reg_moderators, x = moderator, what = "prediction", data = subset(chl_dat, frame_loss == f & opp_name == opp), draw = F, n = x))
      # Label opponent and framing conditions
      marg$opponent <- opp
      marg$frame <- c("Gain", "Loss")[f+1]
      # Bind to results table
      margins <- rbind(margins, marg)
    }
  }

  # Add moderator name to output
  margins$var <- moderator
  
  return(margins)
}
# Apply function to every moderator
chl_opp_marg_frame <- lapply(as.character(moderators_name_chl$var), het_chl_opp_frame) 
# Bind all results into a single table
chl_opp_marg_frame <- bind_rows(chl_opp_marg_frame, .id = "column_label")
# Merge results with moderator labels
chl_opp_marg_frame <- merge(chl_opp_marg_frame, moderators_name_chl, by = "var", all.x = T)

# Create plot
ggplot()+
  geom_line(data = chl_opp_marg_frame, aes(x = xvals, y = yvals, linetype = frame))+
  geom_ribbon(data = subset(chl_opp_marg_frame, !var %in% c("female", "educ_superior")), aes(x = xvals, y = yvals, ymin = lower, ymax = upper, linetype = frame), alpha = 0.3) +
  geom_point(data = subset(chl_opp_marg_frame, var %in% c("female", "educ_superior")), aes(x = xvals, y = yvals), size = 0.2)+
  geom_errorbar(data = subset(chl_opp_marg_frame, var %in% c("female", "educ_superior")), aes(x = xvals, ymin = lower, ymax = upper), width = 0.07)+
  xlab("Moderator")+
  ylab("Predicted prob. of choosing risky option")+
  labs(fill = "", color = "", linetype = "")+
  theme(legend.position = "bottom")+
  facet_grid(moderator~opponent)


### Figure G3 ###
# Create function to generate predicted probabilities as function of moderator
het_chl_hist_frame <- function(moderator){
  # Estimate logit regression
  reg_moderators <- glm(as.formula(paste(
    # Outcome: risk acceptance
    "gamble ~ ", 
    # Moderator fully interacted with treatment variables       
    paste0(moderator, "*frame_loss*historical*hi_val"),
    # Demographic controls excluding ideology to avoid post-treatment bias
    "+ age + female + race_white + educ_superior + santiago")), family = "binomial", data = chl_dat)
  
  ## Estimate predicted probabilities by framing treatment as function of moderator, separately for each opponent condition
  # Estimate at 0 and 1 if moderator is binary; estimate at 10 points if moderator is continuous
  if(moderator %in% c("female", "educ_superior")){
    x <- 2
  }else{
    x <- 10
  }
  margins <- c()
  for(f in 0:1){
    # Estimate predicted probabilities
    marg <- data.frame(cplot(reg_moderators, x = moderator, what = "prediction", data = subset(chl_dat, frame_loss == f), draw = F, n = x))
    # Label framing condition
    marg$frame <- c("Gain", "Loss")[f+1]
    # Bind to results table
    margins <- rbind(margins, marg)
  }
  
  # Add moderator name to output
  margins$var <- moderator
  return(margins)
}

# Apply function to every moderator
chl_hist_marg_frame <- lapply(as.character(moderators_name_chl$var), het_chl_hist_frame) 
# Bind all results into a single table
chl_hist_marg_frame <- bind_rows(chl_hist_marg_frame, .id = "column_label")
# Merge results with moderator labels
chl_hist_marg_frame <- merge(chl_hist_marg_frame, moderators_name_chl, by = "var", all.x = T)

# Create plot
ggplot()+
  geom_line(data = chl_hist_marg_frame, aes(x = xvals, y = yvals, linetype = frame))+
  geom_ribbon(data = chl_hist_marg_frame %>% filter(!var %in% c("female", "educ_superior")), aes(x = xvals, y = yvals, ymin = lower, ymax = upper, linetype = frame), alpha = 0.3) +
  geom_point(data = chl_hist_marg_frame %>% filter(var %in% c("female", "educ_superior")), aes(x = xvals, y = yvals), size = 0.2)+
  geom_errorbar(data = chl_hist_marg_frame %>% filter(var %in% c("female", "educ_superior")), aes(x = xvals, ymin = lower, ymax = upper), width = 0.07)+
  xlab("Moderator")+
  ylab("Predicted prob. of choosing risky option")+
  labs(fill = "", color = "", linetype = "")+
  theme(legend.position = "bottom")+
  facet_wrap(~moderator, ncol = 2)
